home *** CD-ROM | disk | FTP | other *** search
/ Info-Mac 3 / Info_Mac_1994-01.iso / Science / MathPad 2.1.5 / examples / minimize < prev    next >
Text File  |  1993-06-21  |  3KB  |  105 lines

  1. -- Perform multidimensional function minimization.
  2. ~ The user must supply:
  3.  
  4. "f(parms)"  The function to be minimized where
  5.             "parms" is an array of parameters.
  6. "guess"     An array  of initial parameter values.
  7. "deltas"    Initial step size for the parameters.
  8.             If the parameters have different
  9.             scale sizes, "deltas" should define an
  10.             array with a step size for each,
  11.             otherwise a scalar can be used.
  12. "tol"       The termination criteria based on the
  13.             amount of change in parameter values.
  14.             A tol of 1.0e-3 gives about 2 digit
  15.             accuracy.
  16.  
  17. The routine "minimize" steps until a minimum is found.
  18. After minimize is run the following globals are set:
  19.  
  20. "minp"  The solution parameter array
  21. "p"     The final simplex vertices
  22. "y"     The function values at each vertex
  23.  
  24. ~
  25. ------------ Example --------------------
  26. -- fit a curve to data by minimizing the sum of the square of the differences between data and fit.
  27.  
  28. fit(a,x) = a[1]+a[2]*cos(a[3]*x)
  29.  
  30. data=read(xydata); xx[i]=data[i,1]; yy[i]=data[i,2] 
  31.  
  32. err(a) = sum((fit(a,xx[i])-yy[i])^2,i,1,count(data))
  33.  
  34. -- set up globals for minimize
  35. f(a) = err(a)
  36. guess = {80,-80,1}
  37. deltas = .5
  38. tol = 1e-3
  39.  
  40. plot data
  41.  
  42. minimize:;    minp:{81.097,-74.310,0.889}
  43.  
  44. plot fit(minp,X)
  45.  
  46. ----------- minimization routines -------
  47. -- These routines use the "downhill simplex method". A simplex is a N dimensional geometrical figure with N+1 vertices. The algorithim evaluates the function at each vertex and steps the simplex toward the minimum by: 1) reflecting a vertex away from the function's highest point. 2) reflection and expansion away from the highest point. 3) contraction away from the highest point. 4) contraction in all dimensions toward the lowest point.
  48.  
  49. -- step until parameters change by less than tol
  50. minimize = init,
  51.            step,
  52.            step while pchange > tol,
  53.            minp:=psum/m
  54.  
  55. pchange = sum((abs(p[lowest]-p[highest])/
  56.               (abs(p[lowest])+abs(p[highest])))[i],i,
  57.               1,n)
  58.  
  59. -- move simplex one step
  60. step =  rank,
  61.  try(-alpha),                -- reflect from highest
  62.  try(gamma) when ytry ≤ y[lowest], -- keep going
  63.  (ysave:=y[higest],
  64.   try(beta),                 -- contract from highest
  65.   shrink when ytry ≥ ysave) when ytry ≥ y[high]
  66.  
  67. -- find lowest, highest and next highest vertex
  68. rank = lowest:=1,high:=1,highest:=2,i:=1,
  69.    (lowest:=i when y[i]<y[lowest],
  70.     (high:=highest,highest:=i) when y[i]>y[highest],
  71.     high:=i when y[i]>y[high] and i ≠ highest,
  72.     i:=i+1) while i≤m
  73.  
  74. -- move highest vertex through face by fac
  75. try(fac) = ptry:=psum*(1-fac)/n-
  76.                   p[highest]*((1-fac)/n-fac),
  77.     ytry:=f(ptry),
  78.     (y:=replace(y,highest,ytry),
  79.      psum:=psum+ptry-p[highest],
  80.      p:=replace(p,highest,ptry)) when ytry<y[highest]
  81.  
  82. -- shrink entire simplex
  83. shrinkp(k)[j] = .5*(p[k,j]+p[lowest,j]) dim[n]
  84. shrink = kk:=1,      
  85.   ((psum:=shrinkp(kk),
  86.     p:=replace(p,kk,psum),
  87.     kk:=kk+1) when kk≠lowest) while kk≤m,
  88.   y:=fp
  89.  
  90. -- define finite arrays so they can be assigned
  91. sump[j] = sum(p[i,j],i,1,m) dim[n]
  92. initp[i,j]= guess[j]+(deltas[j] when i=j,0) dim[m,n]
  93. fp[i]=f(p[i]) dim[m]
  94.  
  95. init = n:=count(guess), m:=n+1,
  96.        p:=initp,
  97.        psum:=sump,
  98.        y:=fp
  99.  
  100. replace(array,index,newval)[i] = newval when i=index,
  101.                                  array[i]
  102.  
  103. -- tuneable constants
  104. alpha=1; beta=.5; gamma=2
  105.